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Abstract 

This paper analyzes the Krylov convergence rate of a Helmholtz prob- 
lem preconditioned with Multigrid. The multigrid method is applied to 
the Helmholtz problem formulated on a complex contour and uses GM- 
RES as a smoother substitute at each level. A one-dimensional model 
is analyzed both in a continuous and discrete way. It is shown that the 
Krylov convergence rate of the continuous problem is independent of the 
wave number. The discrete problem, however, can deviate significantly 
from this bound due to a pitchfork in the spectrum. It is further shown 
in numerical experiments that the convergence rate of the Krylov method 
approaches the continuous bound as the grid distance h gets small. 



1 Introduction 

The Helmholtz equation plays a central role in seismic imaging, electromagnetic 
scattering and many other applications. For a; e f2 C K'' the equation reads 

(x) = [-A - fc (x)^j u (x) = /(x) , (1.1) 

where the wave number A;(x) depends on the coordinates x and can model, for 
example, the change in refractive index of the material through which electro- 
magnetic waves arc propagating, /(x) models the source of the waves and —A 
is the d-dimensional negative Laplacian. 

In theory many applications need to be solved on an infinite domain, yet 
in practice a numerical solution method must truncate the domain in some 
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way. Therefore, on the finite computational domain the equation is solved with 
outgoing wave conditions on the artificially introduced boundaries. Over the 
years many good outgoing wave boundary conditions have been proposed such 
as Perfectly Matched Layers (PML) [TJ This leads to a spectrum of the 
operator where each eigenvalue has an imaginary part to represent the damping 
of the outgoing waves on the exterior layers. In a similar way in the physics and 
chemistry literature absorbing boundary conditions based on Exterior Complex 
Scaling (ECS) [3l|4j[5] are used for example in break-up problems [6j. 

After discretization the Helmholtz problem becomes a linear system HhUh = 
fh- Due to the negative shift, with a magnitude that depends on the wave num- 
ber k, the matrix A is indefinite. Indeed, the wave number shifts the spectrum 
of — A, which is positive definite, to the left. The eigenvalues corresponding to 
the smooth modes can get close to zero or have a negative real part. These spec- 
tral properties lead to a large condition number and iterative methods perform 
poorly. Efficient preconditioners for the negative Laplacian, such as a multigrid 
method, fail when applied to the Helmholtz problem. A recent review of the 
difficulties of iterative methods for the Helmholtz problem is given in [7]. 

An important step in the improvement of the iterative solution of the Helmholtz 
problem was taken by Bayliss, Goldstein and Turkel [HI [9] with the shifted Lapla- 
cian preconditioner. Instead of approximately inverting the original Helmholtz 
operator with e.g. ILU or a few multigrid cycles as a preconditioning step, the 
Laplacian or positively shifted Laplacian is used as a preconditioner. This pos- 
itive definite operator serves as an approximation of the Helmholtz operator 
and can be efficiently inverted with standard iterative methods. A significant 
extension of this idea was made by the introduction of the complex shifted Lapla- 
cian, a related Helmholtz problem with a complex-valued wave number, which 
makes a better preconditioner and can still be efficiently solved with multigrid 
|101 . The complex- valued wave number prevents that eigenvalues of the pre- 
conditioner come close to zero at any level of the multigrid hierarchy. This is 
particularly useful in the coarse grid correction where diverging numerical res- 
onances can appear when an eigenvalue of a coarse level approaches the origin 

In |13| it was shown that scaling the wave number with a complex value 
has the same effect as scaling the grid distance with a complex value. As a 
result the Helmholtz problem can be efficiently preconditioned by a Helmholtz 
operator discretized on a complex-valued grid. This might be of interest for 
problems where complex-valued grid distances are already used to implement 
the absorbing boundary conditions. This is the case for ECS [6j or PML 

The introduction of complex wave numbers (or grids) avoids the appearance 
of resonances, however, it does not prevent traditional smoothers like oj-Jacobi 
or Gauss-Seidel to be unstable for the smooth modes. In this paper we are in- 
terested in developing a matrix-free method, though we mention that also ILU 
smoothers can be unstable for similar reasons. In [T3] we analyze GMRES as 
a replacement of the traditional smoothers when multigrid is applied to a pre- 
conditioning operator based on complex-valued grids. Numerical experiments 
show that only a few GMRES iterations are needed at every level, as opposed 
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to the results in [T2] where multigrid was used to invert the original Helmholtz 
operator which requires more GMRES iterations at some intermediate levels. 

Note that for the complex shifted Laplacian preconditioner the complex shift 
has a parameter. The choice of the parameter is analyzed in [15] and |16) . 

We mention other promising preconditioning techniques such as Moving Per- 
fectly Matched Layers [T7] , a transformation that turns the Helmholtz problem 
into a reaction-advection-diffusion problem [TH], application of separation-of- 
variables [19] , algebraic multilevel methods [20] , the wave-ray method [STJ [22] , 
and combined complex shifted Laplacian and deflation |23j . 

In this paper we focus the analysis on the wave number dependency of the 
convergence behavior of a preconditioned Krylov subspace method. The pre- 
conditioning operator is the Helmholtz operator discretized on a complex- valued 
grid and it is inverted with multigrid using GMRES as a smoother substitute as 
suggested in [12] and [14] . The paper starts with a review of a one-dimensional 
continuous model problem in Section [2] For this problem the eigenvalues of the 
preconditioned operator are derived analytically and we find that the Krylov 
convergence rate should be independent of the wave number in Section [4] The 
discrete problem, discussed in Section[5j however, does not have this bound. We 
explain the origin of this deviation and give estimates for the different regions in 
the convergence as a function of the wave number k. In Section [6] we illustrate 
the theory with numerical examples. 



2 Model problem 

In this section we formulate a one-dimensional Helmholtz model problem that is 
representative for higher dimensional problems that arise in many applications. 
It is a Helmholtz problem with a constant wave number k on the domain D, = 

[0, 1] e M, 



Hu{x) = 
m(0) = 0; 

u(l) = outgoing wave. 



u{x) =/(x), Vxe(0,l); 



(2.1) 



with a zero Dirichlet boundary condition on the left boundary x — Q and an 
outgoing wave boundary condition on the right boundary x = \. The right hand 
side f{x) represents a localized source term. 



2.1 The Helmholtz problem with ECS 



The outgoing wave boundary condition in (2.1 1 is implemented with exterior 



complex scaling (ECS) [B], an equivalent formulation of the PML technique by 
Berenger [IJ. Therefore the domain is extended to U F = [0, 1] U (1, i?] S K 
after which a complex coordinate transformation is defined as. 



z{x) 



X, 

l + {x 



l)e*' 



[0,1]; 

X e (l,i?]. 



(2.2) 
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We write Rz = z{R) E C for tfie new complex right boundary. This results in 
the domain flUTz — [0, 1] U {l,Rz] G C that is the union of the original real 
domain [0, 1] and a complex line connecting the point 1 to Rz, see Figure [T] In 
this paper we use linear complex scaling by simply rotating the absorbing layer 
over an angle 9^ in the complex plane, but smoother transitions, with a non- 
constant angle, are also possible. Posing a zero Dirichlet boundary condition in 
Rz implies an outgoing wave in the original right boundary x = 1 [2]. 



The Helmholtz problem (|2.1|) translates into 

(2.3) 



Hu{z)= -£,-k^ u{z) =f{z), Vze (0,l]U(l,i?, 



[u{0) = u{Rz) = 0, 

with homogeneous Dirichlet boundary conditions at z{0) — and z(R) — Rz 
Note that the source term f{z) was assumed to vanish outside [0, 1]. 
We define the ECS grid on the complex stretched domain (|2.2|), 



+ {n + l<j<n + m), ^''^^ 

that consists of n intervals of the grid distance h = 1/n followed by m intervals 
of the complex grid distance = {R — \)e'^^'< /m for the complex contour as 
illustrated in Figure [T] We discretize the second derivative operator on the grid 



(2.4) with the Shortlcy-Weller finite difference scheme for non- uniform grids 



, 2 /I /I 1 \ 1 



in grid point j, where hj^i and hj are the left and right grid distance respec- 
tively, and may belong either to the h category or to the hj category. The result 
is a linear system of equations 

HhUh = {-Lh - k^Ih)uh = A, (2.5) 
with a unique solution Uh that approximates the continuous solution u of the 



Helmholtz ec^uation (2.3). The higher dimensional Laplacian A is then con- 
structed with Kronecker products of this one-dimensional discrete Laplacian 
matrix Lh- 

2.2 Spectrum of the discretization matrix 



For the one-dimensional model problem the solution Uh of (2.5) is easily found 
with an exact inversion of the tridiagonal matrix Hh in Equation (2.5). As 
the bandwidth of the matrix grows with the dimension of the problem, so does 
the computational cost of direct methods. Iterative methods need to be used 
instead, such as multigrid and Krylov subspace solvers. The one-dimensional 
model has been analyzed in [13j in order to help in configuring these methods 
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Figure 1: The domain of the model problem (2.3 1 on the continuous domain 
(top) and the discrete grid (bottom). 



efficiently. More specifically their performance depends on the position of the 

T 



eigenvalues of the matrix Hh in the complex plane. Define 7 = tti then the 



eigenvalues of —Lh are the solutions of 

^ tan(2np(^)) cos^ ^ 
^ ^ tan(2mg(t)) cos(g(t)) ^ ^ 

with p{t) = \ arccos(l — |/i^), q{t) = \ arccos(l — |7^/i^). 
Figure [2] shows that the spectrum (•) of ~Lh has a typical pitchfork shape. It is 
bounded in the complex plane by a triangle ioii^2 described by the points = 0, 
ti = 4//i^ and t2 = Starting in the origin = we find eigenvalues along 

the complex line pe^^*^° with p > 0, where 0q is the argument of i?^. It was 
shown in [T^ that these eigenvalues approximate the smallest eigenvalues (x) 
of the continuous Laplacian operator —A that will be derived in Section [3j 
They correspond to the smoothest modes spread over the entire ECS domain 
and we will therefore call them the smooth eigenvalues. At a certain point ti, 
the line of smooth eigenvalues splits up into two branches. One pronounced 
complex branch consists of eigenvalues with associated eigenvectors that have 
their largest components at indices n < j < n + m. Since these eigenvectors 
have nearly-zero components at indices 1 < j < n — 1 that correspond to the 



interior real region of the grid in (2.4), we say that they are mainly located on 
the complex contour of the domain = [Iji?^]. Whereas the other branch 
of eigenvalues in the spectrum lies closer to the real axis and corresponds to 
eigenvectors with their largest components at indices 1 < j < n ~ 1, in other 
words, they are located on the real interior domain O = [0, 1], see also Figure [sj 
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Together with the hne of smooth eigenvalues the latter branch causes potential 
numerical problems as they lie close to the real axis around the points to = 



and ti — 4//i . For the entire Helmholtz operator in (2.5) with a constant 



wave number k the pitchfork shaped spectrum, and the bounding triangle, is 
shifted in the negative real direction over a distance k^. 






real(?L) 



Figure 2: The eigenvalues of the Laplacian discretized on the ECS domain (•) 
lie along a pitchfork shaped figure in the lower half of the complex plane, close 
to the eigenvalues of the Laplacian restricted to the interior real domain (<) and 
the complex contour (>) respectively. The smallest eigenvalues approximate the 
eigenvalues of the continuous Laplacian (x ), until they split up in a point tb (O), 
into two branches with limiting points ti = A/h? and t2 ~ 4/(ft,-y)^. Eigenvalues 
accumulate near the real axis around and ti. 



As the higher dimensional Helmholtz problems are constructed with Kro- 
necker products, these results on the spectrum of the discretization matrix Hh 
are easily extended. Every eigenvalue A of the c?-dimensional Laplacian is a sum 
of eigenvalues A*-'-* of the one-dimensional cases, A = Eti A^'^. This allows us to 
stick the discussion to the basic academic one-dimensional model problem. Note 
that real applications may require more carefully engineered domains with e.g. 
smoother complex stretching, higher order discretization methods or an absorb- 
ing ECS layer on both sides of the domain. These generalizations might have 
an effect on the eigenvalues of the discretization matrix, but the main topology 
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Figure 3: Real part (x) and imaginary part (o) of the three eigenvectors Vo,vi 
and V2 associated to the extreme eigenvalues Aq ~ ioj -^i ~ ti and A2 ~ ^2 of the 
pitchfork shaped spectrum in Figure [2j Eigenvector vq (top) is the smoothest 
and is stretched over the entire domain, interior and exterior ECS contour. 
Eigenvector vi (middle) and V2 (bottom) are highly oscillatory and mainly lo- 
calized on the interior and exterior region respectively of the ECS domain. 
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remains a bounded pitchfork shaped spectrum with the smoothest eigenvalues 
aUgned, close to the continuous case. 

If we assume that the spectrum of the Hclmholtz discretization matrix lies 
inside the triangle totit2 — fc^ in the complex plane, it is straightforward to see 
the main issues for iterative methods. First of all the size of the triangle grows 
as which can be expected with the Laplacian operator involved. This bad 
conditioning destroys the efficiency of Krylov subspace methods. It would not 
necessarily be an issue for a multigrid method, however another difficulty is the 
indefiniteness of the matrix. The negative Helmholtz shift — fc^ drives the upper 
branch of the pitchfork closer towards or even past the origin. This makes the 
coarse grid correction in multigrid highly unstable due to a possible numerical 
resonance at a coarser level as was reported in [THl [H] . A common solution is 
a preconditioned Krylov subspace method where another matrix Mh is defined 
such that, 

can be easily solved instead. The preconditioning matrix Mh is chosen such 
that it is efficiently invertible with a fast multigrid method and such that 
the preconditioned matrix M^^Hh is well conditioned, that is, its eigenval- 
ues are clustered around 1 away from the origin. The complex shifted Laplacian 
^CSL _ _^ _ pf^2 ]-^g^g |-,ggj^ successful choice introduced by Erlangga [TU] 
for Sommerfeld radiation conditions. Simply shifting the Laplacian downwards 
into the complex plane fixes the coarse grid correction in multigrid. In [13] this 
idea was used with ECS boundary conditions, together with the introduction of 
the closely related complex stretched grid (CSG) operator M'^^'^, that is con- 
structed by discretizing the original Helmholtz operator — A — fc^ on a different 
complex stretched domain, 

z(x)-l ''^^''^ ^e[0,l]; 



€'■'^1' + {x - l)e'"-' , x£{l,R\ 

This domain is complex scaled over e*^'' in the interior region [0, 1]; the exterior 
complex contour has the same scaling e'^^ , see Figure |4] The spectrum of the 
discretized operator M^^'^ = L'^^'^ — k^Ih is pitchfork shaped as the original 
Helmholtz operator Hh = Lh — k^Ih, but with the troublesome upper branch 
deeper in the complex plane, see Figure [5] Indeed, back scaling the entire pre- 
conditioning domain over the inner angle dp with e~^^f> returns a regular ECS 
domain with a real interior region and an ECS layer with a reduced angle 9^ — Op. 
Discretizing the Laplace operator on this latter ECS domain gives the scaled 
matrix e'^i'^nL^SG _ ^2j^ ^^^^ M^^'^ = L'j^^^ - k'^Ih must have a pitchfork 
shaped spectrum too, though somewhat more narrow and rotated away from the 
real axis over an angle —29 p. Similar to the complex shift /3 in M'~'^^ , the exact 
choice of the interior scaling angle 9(j determines the performance of multigrid 
on the preconditioner versus the overall convergence rate of the preconditioned 
Krylov subspace method. In [T3] GMRES is suggested as a smoother substitute 
in multigrid which permits small angles 9^ for the preconditioner M'~''^'~^. This 
improves the Krylov subspace convergence significantly. The goal of this paper 
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is to have a better understanding of the spectrum of the preconditioned oper- 
ator M^^Hh, where eventually = M^^'^ will be inverted with a multigrid 
method. 



L 








■ — ■ r% 
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Figure 4: The domain of the CSG preconditioner (solid line) (2.7) differs from 
the original ECS domain (dashed line) in the interior region where it is scaled 
into the complex plane by e^^f . The exterior complex contour has the same 
scaling e*^T as the original ECS domain. 
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Figure 5: The spectrum (•) of the Laplacian discretized on the preconditioning 
domain in Figure [4] is pitchfork shaped too, but with the upper branch rotated 
away from the real axis. 
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3 Eigenvalues of the ID Laplacian on the com- 
plex domain 

In this section we discuss the eigenvalues of the Helmholtz problem formulated 
on an ECS domain as in Figure [T] To this aim we first consider the Laplacian 
on a one-dimensional stretched domain 



z{x) = / q{t)dt, 



(3.1) 



in the complex plane. In order to simplify the discussion in this section we 



purposely use the integral representation (3.1 1 for the ECS domain, as opposed 



to the formulation in (2.2). We are interested in eigenmodes of the Laplacian, 



with Dirichlet boundary conditions Uj{0) = and Uj{z{R)) = 0. For the remain- 
der of this discussion we will drop the subscript j on u and A. After applying 
the chain rule, the equation becomes 



I d 1 d 
q{x) dx q{x) dx 



-A 



u{x) — 0, 



with w(0) = and u{R) = 0. 



The domain is described by Equation (3.1) with 



q{x) 



< a; < r, 



(3.2) 



where r = 1 for the model problem in (2.3). We then have a second order ODE 
with constant coefficients on [0, r] and on (r, R\ and the solutions can be written 
as a linear combination of two fundamental solution. We denote with Ui{x) the 
solution on the first interval and ^2(2^) the solution on the second interval. In 
the point r the solutions of the both subdomains need to be matched with the 
conditions 



ui{r) 

lim,^o^(;^"i(r- 



where the jump condition on the derivative expresses that u{z{x)) needs to 
be continuously differentiable along the transformed domain z{x). Solving the 
equation on each subdomain with boundary conditions mi(0) = and U2{R) — 
leads to 

J Ml (a;) = ^ sin(.T-\/A) , < x < r ; 

\u2{x) = Bs\tl{{x ~ R)-f^/\), r<x<R, 
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where A and B are unknown coefficients. The sohitions ui and U2 have to fulfill 
the matching condition in r 

Ui(r) — U2{r) and u'i{r) = — U2(r). 

7 

After setting A — 1, which can be done without loss of generality, and eliminat- 
ing B with the first condition, the second equation becomes 

/ N sm [rVx] cos {{r - R)-i^/\] 
cos ( r\/A ) ^ 4 ^ ^ = 0. 



0, 



sin {{r - R) 7\/A^ 
After some trigoniometry this leads to the condition 

siri {yX{{R - r)^ + r)^ 

thus to the eigenvalues 



,(i?-r)7 + 
for j e No. 

Note that r + {R — r)^ = R^ = z{R) is the end point of the complex contour 
on which wc have solved the ODE. In this point we have enforced Dirichlet 
boundary condition. Therefore these eigenvalues belong to eigenmodes that are 
standing waves on the domain [0,i?2]. The eigenvalues are independent of the 
details of the complex contour. If we would have taken a more complicated 
contour, with e.g. quadratical scaling instead of linear scaling over a constant 



angle 9^ as in (3.2), the eigenvalues would be the same. 

Note that the discrete problem, discussed in Section [2] only approximates 
the first few eigenmodes along this line. Then, at a certain point along the line, 
the spectrum of the discrete operator will bifurcate into two branches as shown 
in Figure [2] This branch point will be discussed in Section 5.2 



The spectrum of the Helmholtz problem with constant wave number k is 
now ^ 

A,(A;)-(^) -k\ (3.3) 



with j e Nq. These are the eigenvalues of the Laplacian shifted over —k 
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4 Eigenvalues of the preconditioned problem on 
the ID domain 

Let us now look at the eigenvalues of the preconditioning operator M'~^^'^ . It is 
defined on a domain described by 

r < X < R. 
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Let us denote the end point of this complex domain as Rz = p{t)dt. The 
eigenvahies of the preconditioning operator are described by j^n'^/R'^, using 
the results of the previous section. So both for the original problem as for the 
preconditioned operator we have that the eigenvalues lie on a straight line in 
the complex plane. 

Let us assume that for each j the eigenvectors for the domains defined by 
p and q are the same. Then the eigenvalues /i of the preconditioned operator 
{M^^'^)~^ H can be approximated by 



11 j = 



pTT^Rl - e 
fn^/Rl^k^' 



This can be rewritten as 



Rlf'K'^/k^-Rl 



which is the evaluation of a linear fractional transformation LF 
Hz w-R^ ^-j^^ points j^vr^/fc^ i 



with i e Nq. 



(4.2) 



^ (L : w H' 
It is a known property in 



complex analysis that LF maps lines to lines or circles in the complex plane. 
In this case we find that the eigenvalues /ij form a circle in the complex plane 
with radius 



Rz — Rz 



2Q{Rz) 

This circle does not include the origin. Note that the radius of the circle is 
independent of the wave number k. 

However, as the next discussion will show, the assumption that the eigen- 
vectors for a given j on the p and q domains are the same is invalid because 



the eigenvalues of {M'~^^^)~^ H are in fact different from /ij in (4.2). Indeed, 
in order to understand the spectrum we have to solve the eigenvalues of the 
operator 



1 



1 



1 d 1 



p{x) dxp{x) dx J \ q{x) dx q{x) dx 
which is a generalized eigenvalue problem 



1 



1 d 



q{x) dx q{x) dx 
that becomes, after reordering 
I d I d 



u{x) 



1 



p{x) dx p{x) dx 



fc2 u{x) = iiu{x), (4.3) 



u{x) 



1 



1 



q{x) dx q{x) dx 
1 d I d 
q{x) dx q{x) dx 



p{x) dx p{x) dx 



(1 - ^)fc2 



u{x) = 0, 



(l-M)fc' 



u{x) — 0, 



with l/q{x) = \Jl/q{x)'^ — fi/p{xy 
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Proposition 4.1. For the model problems with domains defined by q{r) in (3.2 I 
and p{r) in (4.1), the eigenvalues of the operator in (4.3) are 



with 



1 ( JTT 



7(i?-r) 



Proof. The piecewise constant p and q again lead to a second order ODE with 
constant coefficients for 111(0;) on the interval [0, r] and for 1*2 (s) on the interval 
[r,R]. 

('-(l-/xJj)^Ui-(l-/i)Pu^ ^0 0<a;<r, 
t^)^£2U2- {I- n)k^U2 = r<x<R, 
with iti(O) = and 1*2 (i?) = 0. The solutions ui and U2 are 

{ui{x) ~ Asm ^fc/3y^^j^^a;^ < a; < r, 
U2ix) = Bsm{kj{x- R)) r < x < R, 

that need to be matched by the conditions 

= U2{r); 

= lime^o^(;q:^U2(^ + e)- 



Without loss of generality we can choose A — 1. Requiring continuity of the 
solution in r leads to 



B 



kp./^^Ji)/W^)r 
sin {k^{r — R)) 



Inserting this in the matching condition for the derivatives leads, after some 
trigoniometry, to 



sin kp 



The eigenvalues are the solutions of 



1 — /i 

-r + k"f [R — r) 



0. 



kp 



1-M 



/32 



r + k^ {R — r) — jn 



and we find that 



where 



1' 



1 / JTT 



r \ k 



liR-r) 



□ 
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Corollary 4.2. The eigenvalues of the preconditioned operator (M*^"^*^) lie 
on a parametric curve t : [—^{rj), oo) — >■ C that maps t to 

(t - - 1 



(t-lS(77))2//32- 1' 



with rj = j{R/r — 1). When 3(7) = 0, the curve lies on a circle through 0, P'^ 
and 1. 



Proof. Splitting Sj in Proposition 4.1 into a real and imaginary part leads to 



1 ^^-5R(7)(i?-r)) -zi3(7)(i?-r). 

Define rj — j{R/r — 1), then for each j e Nq there is a i G [—3^(77), +00) C M 
such that Sj — t ^ i{R/r — 1)5(7). In the limit when there is no exterior 
complex scaling, i.e. ^(7) = 0, we can reparametrize hy t = t'^ R and so the 
curve reduces to a linear fractional transformation. It follows that the real line 
is mapped to a circle through the points 0, and 1. □ 

It is important to note that changing the wave number k does not alter the 
parametric curve. Indeed, changing k only modifies the real part of Sj which 
leads to a different particular choice t that gives the position of the eigenvalue 
fij on the curve. This means that there is an upper bound for the condition 
number k = of the preconditioned problem that is independent of k 

and suggests a fast convergence of the preconditioned Krylov subspace method. 

The spread of the eigenvalues on the parametric curve can change as a func- 
tion of k. First note that in the limit j — >■ cx) the eigenvalues /Xj go to /3'^. In 
a similar way, the eigenvalues will accumulate near /J^ as fc — 0. In the other 
case, as k gets larger, the smoothest eigenvalues fij with j <^ k will go to the 
other end of the curve, 

7^(i?/r- 1)2-1 ^ 
j/™o^-' " j^R/r-l)yp^-l^ ' 

leading to a spectrum that is completely spread over the curve. 

We clearly observe this behavior in Figure |6] where the spectrum of the 
preconditioned system is plotted for different values of k for the one-dimensional 



model problem (2.3) with r = 1, R = 1.25, outer ECS angle 9, 



1 



angle for the preconditioner Op — 0.18 w jj. The circles mark the first 80 



eigenvalues of the continuous problem, as given by Proposition 4.1 They lie 
on the parametric curve derived in Corollary |4.2| which is visualized with a 
solid line. For a small wave number k = 0.4, in the upper left subfigure, all 
eigenvalues fij, with 1 < j < 80 lie close to (3^. In the upper right plot with 
k = 6.4 we see that /ii and /i2 almost reach the other end of the curve, while 
the remaining eigenvalues fj,j with j > 3 still lie closer to /J^. Finally in the 
lower two subfigures the eigenvalues are more spread along the curve as k grows 
larger. 
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In a similar way the 2D eigenvalues, or for any higher dimension, will be 
bounded by a parametric curve that is independent of the wave number. In 
contrast to the one-dimensional case the 2D spectrum will fill up the region 
bounded by the curve with eigenvalues. 



k = 0.4 k = 6.4 




0.4 0.6 0.8 1 0.4 0.6 0.8 1 

real(|i) real(f^} 



Figure 6: First 80 Eigenvalues of the one-dimensional preconditioned system 
for different values of the wave number k. For the continuous operator (o) the 
spectrum lies on a parametric curve (red line) in the complex plane. When k 
increases the eigenvalues spread over the curve. The spectrum of the discrete 
operator (•) approximates the continuous case for k = 0.4, 6.4 and 16.4 leading 
to an initally bounded condition number in Figure |8] For k = 26.4 the discrete 
eigenvalues lie outside the curve resulting in a larger condition number. For 
higher dimensional problems the region bounded by the curve will fill up with 
eigenvalues. 



5 Discrete operator 

5.1 Deviations from the continuous problem 

However, when the problem is discretized, with for example finite differences, 
the Krylov convergence rate can differ significantly from the bounds predicted 



15 



by the analysis of the continuous problem in Section [3] and [4] 

Indeed, Figure [7] shows the number of GMRES iterations to solve the Helmholtz 
problem with discretization matrix Hh, preconditioned with the complex shifted 
grid matrix M^^'^ that is exactly inverted, as a function of the wave number 
k. The results are for a ID (solid line) and a 2D (dashed line) problem with 80 
grid points per dimension, n = 64 interior grid points and to = 16 additional 
points for the ECS layer. For the 2D problem we recognize three regions in the 
convergence rate. First, between fc = and k sa 16.4, the number of iterations is 
ramped up from a few to about 18. The second region is between k w 16.4 and 
k « 21, where the number of iterations remains constant. Finally from k — 21 
on the number of iteration rises again. This is in contrast to the analysis of the 
continuous problem that predicts a convergence rate that is independent of /c, 
once it is large enough. These observations can be explained with the help of the 
spectrum of the discrete preconditioned system Mf^^Hh of the one-dimensional 
problem shown in Figure [6j 

For the wave numbers in the first region, < fc < 16.4, the spectrum of the 
discrete problem (•) is a good approximation for the spectrum of the contin- 
uous continuous problem (o), except for one or two spurious eigenvalues. As 
a consequence, the spectrum lies along the parametric curve (solid line). The 
initial growth in the convergence behavior is due to the fact that the eigen- 
values are accumulated near for small k and start spreading over the curve 
towards 1 for larger k, leading to an increasing condition number k — '"'^'^j 
where 1 < j < 80. In Figure [8] the condition number is plotted as a function 
of the wave number fc. Indeed, first the condition number is close to 1 until 
eigenvalue fi^o moves significantly along the curve reaching the point with the 
smallest possible absolute value for fc « 2.5. The next eigenvalue /iyg gives rise 
to a second peak around fc « 5 when it reaches that same minimal point on 
the curve. As the spectrum starts spreading more equally over the curve, see 
e.g. lower left subfigure in Figure [6j the condition number seems to converge 
to /« « 2.5 and the number of GMRES iterations stagnates around 18 for the 
second region fc 16.4 to fc w 21 on Figure [7] 

However, from fc w 16.4 on the discrete condition number starts behaving 
differently from what the continuous eigenvalues predict. Around fc w 21 it even 
grows above the upper bound given by the curve. This is because the eigenval- 
ues of the discrete preconditioned system start deviating from the continuous 
spectrum, see e.g. lower right subfigure in Figure [6j Whereas the continuous 
eigenvalues remain on the curve, the discrete spectrum grows outside the curve 
from a certain critical wave number. This divergence between the continuous 
and the discrete preconditioned spectrum finds its origin in the spectrum of the 
original Helmholtz operator. At the end of Section |3] we mentioned that only 
the first few eigenvalues Xj in (3.3) of the continuous Helmholtz operator H are 
well approximated by the eigenvalues of the discretization matrix Hfi given by 
(2.6). In Figure [2] we see that the continuous spectrum (x) lies along a line in 
the complex plane, whereas the discrete spectrum (•) branches from this line 
at a certain point € C. The location of this point will be an indicator for the 
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start of the third region in the convergence behavior of GMRES. For the 2D 
problem the three different convergence regions are more pronounced than for 
the ID problem because the region bounded by the curve gets filled with extra 
eigenvalues. 
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Figure 7: Number of GMRES iterations to solve the one-dimcnsional (solid line) 
and two-dimensional (dashed line) Helmholtz problem as a function of the con- 
stant wave number k with 80 grid points per dimension. For the 2D problem 
three regions are clearly distinguishable. Initally the number of iterations in- 
creases with k until a stagnation of 18 iterations is reached around k — 16.4. 
For wave numbers larger than fc w 21 the number of iterations starts increasing 
again. Four values of k for which the spectrum of the preconditioned system is 
shown in Figure are marked (□). 



5.2 Predicting the branch point 

Next we derive an explicit formula for the approximate position of the branch 
point in the spectrum of the discretization matrix, see Figure [2j From this point 
on, the eigenvalues of the discrete problem start to deviate significantly from 
the eigenvalues of the continuous problem. This branch point will predict from 
which k on we can expect a rising cost of the numerical solution method. A 
surprising result of this section is that this branch point does not shift with the 
order of the discretization. 
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Figure 8: Condition number of the continuous (o) and discrete (•) precondi- 
tioned system. The condition number of the continuous operator converges to a 
fixed value determined by the parametric curve in Corohary |4 . 2| whereas the dis- 
crete operator shows diverging behavior from a certain wave number k w 16.4. 
Four values of k for which the spectrum of the preconditioned system is shown 
in Figure |6] are marked (□). 



Proposition 5.1. Consider a one- dimensional grid as in (2.4) defined on an 



ECS domain consisting of two parts, [0, r] for the interior with n grid points, 
and the complex interval [r, Rz] = [r, z(R)] for the exterior part with m grid 
points. The smallest eigenvalues of the negative Laplacian discretized on this 
ECS grid —L^, with zero Dirichlet conditions at the boundaries and Rz, lie 
along the complex line 



t{p) = 



with p > 0, 



close to the eigenvalues of the continuous operator tj — (^j^^ with j G Nq. 
For larger eigenvalues the spectrum of —Lh splits into two branches around the 
point tb — (^'R'^ with 



Pb 



R-r 



Rz — R 
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where = z{R) and W{.) is the Lambert-W function. 



Proof. The spectrum of the Laplacian discretized on an ECS grid has a pitchfork 
shape. In order to find the point where the pitchfork spHts we start from 



the condition (2.6) that has the eigenvalues of the discrete Laplacian —L^ as 



solutions. In the rest of this discussion we will consider the scaled Laplacian 



L = —h Lh. The eigenvalues of L are again the solutions of (2.6 1, but with 
p{t) = \ arccos(l — |), q{t) = \ arccos(l — |7^) instead. The condition (2.6) is 
equivalent to 

Fi{t) = sin(2np(t)) cos(2mg(t)) cos(g(t)) + cos(p(t)) cos{2np{t)) sm{2mq{t)) = 0, 

(5.1) 

Since we are interested in the smallest eigenvalues of L we can use the ap- 
proximate condition 



F, 



(t) = sin ^(n + m^) Vt^ — - sin (^nVt^ cos (^mjVt^ tan ^"^^ ~ 



with £ = 7 — 1. 

This is easily derived from (5.1) by using the Taylor series p{t) — ^ 



C(|t|3/2) ^j^j ^(-^^ = 7^ + 0{\t\^/^), for \t\ < 1 and substituting 7 = 1 + e. 
Moreover, since |e| < 1 for realistic exterior complex scaling with an ECS angle 
< f , we have used cos (7^) = cos (^^^ - i sin (^^) Vie + 0{\te^\). 

We will now look at the evaluation of function F2 along the complex line 

t(p) = ( ^ with p > 0. 

\n + J 

It returns real numbers between —1 and 1 for the first term of ^2, and complex 
numbers for the second term. The latter term is small, for small p, and thus 

the roots of F2 will approximately be the roots tj — (^ „^^^ j , with j G No, of 

the first term. The eigenvalues will branch from the line t{p) when the second 
term of F2 becomes more important. This is when, 

A ■ f P \ f P \ f P \ P , . 

- sm n cos m-f tan — r e\ ~ 1, 

2 \ n + rri'j J \ n + m"f J \2(n + m^) J n + m"f 
-l|Bin(,(l-4^))(-^)%i, 
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and after using the identity Rz = r + {R — r)7, 



<4>|- Sm p I 1 



Rz 



ph 

R~z 



h\V~e\ 



" it 



where W{c) is the Lambert-W function, defined such that c = W{c)e^'^''\ and 
evaluated in 



hWe\ 



1 R-r 



Rz y Rz — R 



with £ = 7 — 1 = . The point t;, along the line where the pitchfork 

splits into two branches is now approximately given by 



Pb 



R-r 



1 



Rz y Rz — R 



Pb 



n + 



Pbh 

Rz 



So for the eigenvalues of the unsealed operator L we have tb = 



□ 



The point ti, = 



predicts the point in the spectrum of the discrete 



Laplacian where the pitchfork splits. The smallest eigenvalues lie close to 
with j £ No such that jn < pi,. This is illustrated in Figure |9] where the 
32 smallest eigenvalues are plotted for three different grid sizes n, together 
with the branch point predictions if,, for the ECS domain [0,r] U {r, Rz] = 
[0, 1] U (1,1 + 0. 256"/^]. 



Figure 10 shows the distance to the origin \ti,\ of the predicted branch point 
as a function of the interior grid size n, the other domain parameters are fixed. 
The branch point was also detected experimentally by measuring the deviation 

from the line t(p) = (^^^ with p > 0. As the grid size n increases, the tail 
of the pitchfork grows proportional to the square of the Lambert W-function, 
\tb\ ~ W{n)'^, and not according to the order of discretization. 

The spectrum of the indefinite Helmholtz discretization matrix Hh is achieved 
by shifting the pitchfork spectrum of the Laplacian to the left in the complex 
plane over a distance determined by the wave number k. We can now pre- 
dict the value k for which the eigenvalues of the discrete preconditioned system 
M^^Hh start growing outside the parametric curve of the continuous precondi- 
tioned spectrum. Since tb marks the point in the pitchfork shaped spectrum of 



20 



-1000 



-2000 



-3000 



-4000 



• n 


= 32 


• n 


= 512 


n 


= 8192 



1000 2000 3000 4000 5000 6000 
real(?i) 

Figure 9: The 32 smallest eigenvalues (•) of the discretized Laplacian for n = 32, 
n = 512 and n ~ 8192. The branch point tb (O) in the spectrum moves further 



in the complex plane as predicted by the formula in Proposition 5.1 
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Figure 10; The absolute value of the branch point tf, of the pitchfork in the 
spectrum of the discretized Laplacian for n — 2^ with j = 5, . . . , 20. As predicted 



by the formula (•) in Proposition 5.1 measured experimentally (o) and with a 
higher order scheme in the turning point r (□). As the grid size n increases, 
the length of the typical line of smooth eigenvalues grows proportional to the 
square of the Lambert W-function, |tf,| ~ W{n)'^, and not as the order of the 
discretization. 
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Hh where the badly-approximated continuous eigenvalues lie on the right hand 
side, we expect this to happen when tt is of the same order of magnitude as 
the smooth eigenvalues, this means when \tb\ = k^- In Figures [t] and [s] the crit- 
ical wave number kf, = -y/ftbl = 17.9327 is marked. Indeed, for wave numbers 
k < kj, the eigenvalues of the continuous preconditioned operator lead to a good 
prediction of the effective discrete condition number and the rise in GMRES 
iterations is a direct consequence of the spreading along the parametric curve. 
After a short region of stagnation where the curve is densely filled with eigen- 
values, the number of iterations increases again, however, this time due to the 
growing deviaton of the discrete eigenvalues from the continuous eigenvalues. 
In the next section we will study the convergence behavior for wave numbers in 
this third region with numerical experiments. 



6 Numerical experiments 

In this section we analyze the performance of a Krylov subspace method applied 
to a Helmholtz problem with constant wave number fc in a square domain Q = 
(0, 1)^ with a centered point source 

ifx = l/2 = y, (g^^^ 
elsewhere. 

This two-dimensional problem is built with Kronecker products of the one- 




dimensional model problem (2.3) with outgoing wave boundary conditions in 
every direction. All boundaries are therefore extended with an ECS layer with 
an angle 9^ = ^ to absorb outgoing waves. The results from the previous 
sections are still useful if we take into account that the spectrum of the two- 
dimensional operator is the set of all possible sums of two eigenvalues of the 
one-dimensional case. This means the spectrum of the discretization matrix Hh 
looks like a sum of pitchforks now, as discussed in Section [5T| with three points 
close to the real axis t — —k'^ , t — — k'^ and t = ^ — fc^. These eigenvalues 
correspond respectively to the smoothest mode on the domain, the eigenmode 
which is oscillatory in one dimension and smooth in the other and finally a mode 
that is oscillatory in both dimensions. 

The complex stretched grid preconditioning matrix M^^'^ is constructed by 
discretizing the problem on a complex stretched grid with a small inner angle 
Op = 0.18 « Yf- This ensures that for every level in the multigrid hierarchy the 
eigenvalues are bounded away from zero. The angle for the outer ECS layers is 
kept dX 0^ = ^ as in the original Helmholtz problem as illustrated in Figure |4j 
For a detailed description of the spectrum of this preconditioning matrix we 
refer to [13] and 

As a consequence the preconditioning matrix Mj^^^ can efficiently be in- 
verted with a multigrid method with either three steps of GMRES as a smoother 
substitute, denoted as GMRES (3), or a specific polynomial smoother as sug- 
gested in [M]. However, because the smoother can differ each application the 
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actual preconditioner is not the same in every outer Krylov step. Therefore 
FGMRES, the flexible GMRES method [H], is used as outer Krylov subspace 
methods. We discuss the performance of preconditioned FGMRES before con- 
vergence to a residual norm of order 10~^. The preconditioning matrix M^^*^ 
is approximately inverted with one V(l,l)-cycle with GMRES(3) as smoother. 
The experiments are all run in Matlab"'"'^ on two quad core Intel"'"'^ Xeon 
CPUs (E5462 @ 2.80GHz). 

Figure [Tl] shows the convergence results for preconditioned FGMRES to 
solve the 2D Helmholtz problem (6.1 ) with a residual norm below 10^^ for wave 
numbers ranging from fc = 15 to fc = 180. For these wave numbers the smoothest 
eigenvalues of the discretization matrices have a negative real part. Each curve 
shows the same experiment for a fixed grid size n in one dimension, the grid size 
of the ECS layer is related as m = n/A. Just like the wave numbers the grid 
sizes are purposely chosen over a wide range as well, from n = 16 to n = 2048 
in one dimension, in order to expose the full effect of the preconditioning on 
the convergence behavior of FGMRES. A discussion on the physical accuracy 
of the grid sizes lies not within the scope of this analysis. As we are merely 
interested in the convergence behavior of Krylov subspace methods we explain 
these curves as a function of the increasing wave number k. 

For each grid size both the convergence rate and the number of iterations 
grow initially as a function of A; up to a peak where k ^ | = 2n. This corre- 
sponds to the first critical point t = — k"^ where the pitchfork in the spectrum 
of Hh nearly touches the real axis. These eigenvalues correspond to eigen- 
modes that oscillate rapidly in one direction while they are smooth in the other. 
There is now an eigenvalue of Hh near the origin and the preconditioned system 
M^^Hh obviously suffers from this too. As the wave number k increases more, 
the pitchfork is shifted further to the left in the complex plane. The convergence 
improves slightly until the second critical point t = ^ — k"^ comes too close to 

the origin for k ss = 2^/2n. After this, the spectrum has completely shifted 
into the negative real part of the complex plane, making the Helmholtz matrix 
negative definite. This is observed on the curves as a sudden improvement in 
convergence. As a reference these experiments were repeated for the smallest 
grid sizes with regular GMRES and an exact inversion of the preconditioner 
^CSG ^ in order to eliminate the effect of the approximate multigrid inverse. In 
Figure [T2| the described convergence behavior is then more pronounced. 

Only for the two smallest grid sizes this behavior lies completely within the 
tested range of wave numbers k. Indeed, in Table [T] the critical wave numbers 
are ki = = 2n and k2 = = 2^/2n are listed for the different grid sizes in 
the experiments. The first critical wave number fci has the worst performance 
because the spectrum of iJ^ reaches its most extreme indefiniteness, with the 
pitchfork perfectly spread over the negative and positive real part of the lower 
half of the complex plane. For larger k the spectrum tends more to negative 
definiteness. The second peak is right before it turns completely negative definite 
in the second critical wave number k2 after which the convergence rate obviously 
drops drastically. 
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n 


16 


32 


64 


128 


256 


512 


1024 


2048 


ki 


32 


64 


128 


256 


512 


1048 


2048 


4096 


k2 


45.3 


90.5 


181.0 


362.0 


724.1 


1448.2 


2896.3 


5792.6 


kb = \/\tb\ 


12.8 


13.8 


20.3 


26.2 


35.7 


40.6 


50.0 


53.1 



Table 1: The critical wave numbers fci and ^2 for the discrete Helniholtz prob- 
lem with constant k for different interior grid sizes n. For these values of k 
the preconditioned Krylov method reaches its worst performance as we see in 
Figures [Tl] and [T2| The absolute value of the branch point indicates a wave 
number kb = \/\tb\ that marks the end of an early plateau in the convergence 
behavior discussed in Section 5.1 



For a realistic physical solution the grid size should be large enough in order 
to represent the wave accurately. Higher wave numbers k require finer meshes 
|25| . This means that in practical circumstances only the region on the curve 
long before the first peak is important, where the spectrum is only slightly 
negative definite. In that region we still profit from the fact that for fc — the 
eigenvalues of the prcconditioner accumulate into a single point (see Section 5.1 ). 
There we see a rise and the stagnation of the number of iterations into a plateau 
with the branch point tb as indicator. However, on Figures 11 and 12 the range 
of wave numbers k is too wide to clearly uncover this initial effect as in Figure [Tj 



7 Discussion and Conclusions 

In this paper we have analyzed the convergence rate of a multigrid precon- 
ditioned Krylov solver for Helmholtz problems with absorbing boundary con- 
ditions. The multigrid method inverts a prcconditioner that is a Helmholtz 
operator discretized on a complex-valued grid rather than on a real grid. This 
prcconditioner is comparable to the complex shifted Laplacian. The multigrid 
method uses GMRES(3) as a smoother at each level. 

To understand the Krylov convergence, we have proposed a model problem 
with a Dirichlet boundary condition on one side and an outgoing wave boundary 
condition at the other. The outgoing boundary condition is implemented with 
exterior complex scaling (ECS) that extends the domain with a complex- valued 
contour. This model problem is representative for the implementation of ab- 
sorbing boundary layers in various applications such as ECS, which is often used 
in chemistry and physics or PMLs, which are frequently used in engineering. 

We have analyzed this model both in a continuous and a discrete way. For 
the continuous problem we have found that the spectrum of the preconditioned 
operator lies on a curve in the complex plane which is bounded away from 
zero. This leads to an expected Krylov convergence rate that is bounded for all 
wave numbers. For small wave numbers the convergence rate is faster since the 
eigenvalues accumulate to a single point. 

In the discrete problem the spectrum behaves similarly to the continuous 
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problem for small wave numbers. However, for larger wave numbers the spec- 
trum can deviate significantly. This finds it origin in the properties of the 
discrete Helmholtz operator that has a pitchfork in the spectrum, where only 
one of the arms of the pitchfork approximates the spectrum of the continuous 
operator. These deviations destroy the nice convergence expectations given by 
the continuous operator. The distance to the origin of the predicted point tb 
where the spectrum bifurcates grows only very slowly as a function of the num- 
ber of grid points. Numerical experiments show that also varying wave numbers 
result in similar pitchfork shaped spectra with branching points that can still 
be fairly well estimated using Proposition |5.1[ 

As a rule of thumb the Krylov convergence is bounded when fc^ is smaller 
than the absolute value of the branch point tb- There we can expect a bounded 
convergence rate. For wave numbers k larger than this branch point the number 
of iterations rises until fc^ « 4//i^, where the number of iterations is maximal. 
We have a second but milder peak at ~ 8//i^. From then on the spectrum is 
negative definite and the convergence rate drops rapidly and only a few iterations 
are required to solve the system. 

We conclude that there is no overall /c-independent convergence rate, yet the 
number of iterations diminishes as the number of grid points is increased. To 
further improve the convergence rate of the iterative method it is possible to 
engineer the parameters of the absorbing layer as a function of the wave number. 
For large wave numbers k we do not need a large ECS grid or a large rotation 
angle to absorb the wave. Adapting the parameters can reduce the number of 
iterations to solve the problem. 

Although the current analysis is for constant wave numbers fc, we believe 
many results will still be valid when the wave number varies over space. Indeed, 
a space-dependent k{x) will only affect the smoothest eigenvalues while the ex- 
treme values that determine the diverging behavior depend on the grid distance 
and remain the same. 

There still remain important challenges in the development of an efficient 
solver for the Helmholtz problem with space-dependent wave numbers based 
on complex stretched domains. In numerical experiments with strongly space- 
dependent wave numbers that allow evanescent waves we have seen serious de- 
teriorations of the convergence rate [H] . This is caused by the multigrid coarse 
grid correction on levels too coarse to resolve the evanascent waves. This is a 
subject of future research. 
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(a) The convergence rate of FGMRES, averaged over all iterations. 




k 

(b) The number of FGMRES iterations. 

Figure 11: Convergence results of FGMRES as a function of the wave number 
k for the homogeneous 2D Helmholtz problem. The preconditioner is approxi- 
mately inverted with one V(l,l)-cycle. The lines represent different interior grid 
sizes 16 X 16, 32 x 32,. . . , 2048 x 2048. 
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(a) The convergence rate of GMRES, averaged over all iterations. 




(b) The number of GMRES iterations. 



Figure 12: Convergence results of GMRES as a function of the wave number 
k for the homogeneous 2D Helmholtz problem. The preconditioner is exactly 
inverted. The lines represent different interior grid sizes 16x 16, 32x32,. . . , 512x 
512. The convergence patterns resemble those in Figure 11 for the FGMRES 
case where the inverse of the preconditioning matrix is approximated with one 
V(l,l)-cycle. The discussed effects are more pronounced. 
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